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Current theoretical treatment of mode splitting and scattering loss resulting from sub- wavelength 
scatterers attached to the surface of high-quality-factor whispering-gallery-mode microresonators 
is not satisfactory. Different models have been proposed for two distinct scatterer regimes, i.e., 
a-few- and many-scatterers. In addition, many experimental results seem difficult to understand 
within the existing theoretical framework. Here we develop a unified approach that applies to 
an arbitrary number of scatterers, which reveals the applicable conditions and the limits of the 
existing theoretical models. Moreover, many new understandings on mode splitting and scattering 
loss have been achieved, which are supported by numerical and experimental evidences. Such a 
unified approach is essential for the fundamental studies as well as the practical applications of 
mode splitting and scattering loss in high-quality-factor whispering-gallery-mode microresonators. 



I. INTRODUCTION 

Whispering-gallery-mode (WGM) microresonators 
have attracted a lot of research interest due to their 
high quality factors (Qs) and microscale mode volumes 
[Hi]. The high Q/V factor, i.e., the so-called Pureed 
factor pQ, enables strong light-matter interactions and 
is the essential reason behind the wide applications of 
WGM microresonators including low-threshold lasing 
[3] , low-power optical modulation 4 , single-nanoparticle 
sensing [5J [B], ultrasensitive micromechanical displace- 
ment detection \7[, as well as the fundamental studies 
on cavity quantum electrodynamics [8 . Because of 
the structural symmetry, ideal WGM microresonators 
usually exhibit degeneracies in their resonant modes. 
For example, clockwise (CW) and counterclockwise 
(CCW) WGM modes are supported in microtoroid and 
microdisk resonators with identical mode properties 
(e.g., resonance frequency and linewidth). Reflected in 
transmission, only one single resonance is observed. 

The degeneracy in the resonant modes can be lifted 
by destroying the structural symmetry, either done in- 
tentionally such as by introducing nanoparticles to the 
surface of microresonators [9[ [10], or caused by imper- 
fect fabrications which distort the structure |llj . In con- 
sequence, the two degenerate modes (i.e., the CW and 
CCW modes) will couple to each other and a doublet ap- 
pears in the transmission, a phenomenon termed as mode 
splitting [12] . The structural defects can also couple 
the confined WGM modes to free-space radiation modes, 
generating scattering loss and thus linewidth broadening 
to the WGM modes [9]. 

Mode splitting and scattering loss have been investi- 
gated in many different works [5HTT1 [TMT5] . Based on 
their applications and the number of scatterers, previ- 
ous works can be categorized into two distinct scenarios. 
In the first scenario, small nanoparticles are introduced 
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to the surface of high-Q microresonators, with focused 
applications such as strong light-matter interactions and 
nanoparticle sensing [5] [SJ [S] . The number of nanopar- 
ticles is usually limited to a few, and reasonably good 
agreements between experimental observations and de- 
veloped models have been achieved. The other scenario 
considers sub-wavelength scatterers that are intrinsic to 
microresonators, such as surface roughness caused by im- 
perfect fabrications which is of fundamental importance 
for a thorough understanding of high-Q microresonators 
[llj . In such cases, the number of scatterers is typically 
on the order of hundreds or even thousands, and a differ- 
ent approach has to be taken to study the mode splitting 
and scattering loss [TT] . 

As we shall detail in section II, existing theoretical 
works developed for the above-mentioned two scenarios 
have only achieved partial success. Roughly speaking, 
they can explain the mode splitting quite well but have 
major problems in predicting the corresponding scatter- 
ing loss even qualitatively. Though the models developed 
for the a-few-scatterer case have claimed possible exten- 
sions to the many-scatterer scenario [US [14], we will show 
that they have neglected some fundamental effects and 
such an extension will not be successful. On the hand, 
the approach developed for the many-scatterer case has 
its own shortcomings. For example, in experiments the 
two split modes usually exhibit different linewidths, a 
fact that cannot be explained by the model shown in 
Ref. |llj in a self-consistent manner. Moreover, many 
other experimental observations, such as the azimuthal- 
order variations of the mode splitting and scattering loss 
within the same radial mode family in an individual mi- 
croresonator |16j . are difficult to understand within that 
framework. 

In this paper, we will develop a unified approach that 
applies to an arbitrary number of scatterers. Our results 
are closely compared to those of the existing approaches 
for the two distinct (i.e., a-few- and many-scatterer) sce- 
narios , and we show conditions under which our model 
can be reduced to the existing models in their respective 
regimes. Moreover, our work has provided new under- 
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standings on the mode splitting and scattering loss in 
high-Q WGM microresonators, which are supported by 
numerical studies and experimental results. For example, 
there is an intuitive belief that in the presence of mode 
splitting, the eigenmode exhibiting a lower resonance fre- 
quency of the doublet would also incur a stronger scat- 
tering loss from positive dielectric perturbations, which 
is true for the single-scatterer case and also consistent 
with existing theoretical models O HI] (see detailed dis- 
cussion in section II). However, we will prove that such 
perception is generally invalid. Another example is that 
our model predicts when the two eigenmodes overlap in 
the resonance frequency (i.e., no mode splitting), gen- 
erally their scattering loss rates are different. There- 
fore, the two eigenmodes are not degenerate in the strict 
sense, a fact that has been overlooked in previous stud- 
ies [TOl HI]. Furthermore, we apply our model to the 
fabrication-induced surface roughness in high-Q WGM 
microresonators, which unveils that the mode splitting 
for different azimuthal orders are statistically indepen- 
dent. Hence, for the same radial mode family in an in- 
dividual microresonator, strong variations of mode split- 
ting are possible. In addition, we have shown that the 
scattering loss of the same radial mode family can exhibit 
more than 30% variations among different azimuthal or- 
ders, while initially one might expect a uniform scattering 
loss rate for these modes. This in fact solves one mys- 
tery that often confuses people working on high-Q mi- 
croresonators, that is, in scattering-loss-limited microres- 
onators, the extremely high intrinsic Q could only be ob- 
served for one azimuthal order while the intrinsic Qs of 
the rest modes could be much lower 



II. 



EXISTING MODELS VERSUS 
EXPERIMENTS 



In this section, we will review major existing models on 
the mode splitting and scattering loss in high-Q WGM 
microresonators. As mentioned, currently there are two 
different approaches working at two distinct regimes of 
scatterers. By comparing the theoretical results pre- 
dicted by these models to experimental (or numerical) 
observations, it is shown that the two approaches are 
only partially successful, in the sense that they either are 
not self-consistent or fail to agree with some of the ex- 
perimental results. Therefore, a unified approach, which 
could provide a full understanding on the mode splitting 
and scattering loss, needs to be developed. 

The first approach considers a coupled system consist- 
ing of the CW and CCW WGM modes as wells as free 
space radiation modes, with their interactions assisted 
by each individual scatterer. The single-scatterer case 
has been well studied in high-Q microtoroid resonators 
[9j [13] . Two standing- wave modes, being symmetric and 
anti-symmetric combinations of the CW and CCW trav- 
elling modes, appear in the resonance spectrum. The 
symmetric mode has a nonzero field overlap with the 



scatterer, resulting in a red shift in its resonance fre- 
quency and a broadening in its linewidth (we assume pos- 
itive dielectric perturbations from scatterers throughout 
the paper unless specified). On the other hand, the anti- 
symmetric mode has a zero field overlap with the scat- 
terer; therefore, its resonance frequency and linewidth 
stay the same as those of the CW (CCW) mode with- 
out the scatterer. These two resonances are illustrated 
in Fig. 1(a), with w+ and 7+ (j-) denote the reso- 
nance frequency and loss rate of the eigenmode that has a 
high (lower) resonance frequency of the two split modes 
(i.e., w+ > UJ-). The loss of a high-Q microresonator 
can come from many sources, such as scattering loss due 
to scatterers and absorption loss from material absorp- 
tion. Since we are primarily interested in scattering loss 
in this paper, other loss mechanisms are not considered 
unless specified. As illustrated in Fig. 1(a), for the single- 
scatterer case, we have 7+ < 7_. The two-scatterer sce- 
nario has also been explored in Ref. [TU], which has ex- 
perimentally demonstrated that the two eigenmodes can 
either have no mode spitting (cj + = lj_) or a symmet- 
ric (7 + = 7_) or an asymmetric (7+ < 7_) doublet, 
depending on the relative position of the two scatter- 
ers. These experimental results can be analyzed by a 
generalized model presented in Ref. [M] , which considers 
multiple scatterers that are well separated apart so that 
the contribution from each scatterer can be considered 
to be independent with each other. However, when the 
scatterers are closely spaced, this independent-scatterer 
approach fails to predict correct results. For example, for 
N identical scatterers, the model in Ref. p3] gives 
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where uj c is the originally degenerate resonance frequency 
of the WGM modes before the introduction of scatterers; 
Go and r are parameters (both positive) characterizing 
the resonance frequency shift and linewidth broadening 
caused by an individual scatterer, respectively; k is the 
wavenumber of the WGM mode; and x n is the projection 
of the nth scatterer's position on the WGM's wave trav- 
elling direction. From Eq. (2), one can infer that for an 
arbitrary number of identical scatterers, 7+ < 7-- This 
conclusion seems valid, by arguing that the resonance 
corresponding to uj_ has a lower resonance frequency be- 
cause of the more field overlap with the scatterers, which 
is also responsible for a stronger scattering loss. How- 
ever, such intuitive belief is not generally true. Here we 
consider one extreme example in Fig. 1(b), where identi- 
cal scatterers have uniformly covered the outer surface of 
a microresonator except for one vacancy. Assuming the 
scatterers and the microresonator share the same dielec- 
tric constant, the resulting structure can be regarded as 
a larger-radius microresonator with a negative-dielectric- 
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FIG. 1. (a) Transmission response of a single dielectric scat- 
terer on the surface of a microresonator. (b) Transmission 
response of numerous dielectric scatterers which uniformly 
cover the surface of a microresonator except for a vacancy. 



constant scatterer at the vacancy point. Consequently, 
the mode that has a nonzero field overlap with the scat- 
terer will incur a blue shift in the resonance frequency 
(relative to the resonance frequency of the larger-radius 
microresonator) as well as a linewidth broadening from 
the scattering. Hence, j + > j_, contrary to the result 
from Eq. (2). Another problem with the independent - 
scatterer approach is that it has to track the position of 
each scatterer, which makes it impractical for problems 
such as surface roughness caused by imperfect fabrica- 
tions in high-Q microresonators, where the positions of 
scatterers are random and only statistical information is 
available. 

The second approach, developed mainly for the sur- 
face roughness present in high-Q microresonators, em- 
ploys an intuitive physics model (which is essentially a 
phenomenological model) to describe the mode splitting 
and uses the volume current method to obtain the scat- 
tering loss Here, we use the microdisk resonator 
as an example to give an introduction to this approach. 
For an isolated microdisk resonator (i.e., no external cou- 
pling), we have [T3] 

dClccM, ( . .. 7CCW\ ■ Q / \ 

■ y-iu c -+- !nw ccw -\ — j a ccw -r- ip 

^ : - (iu c + iAuj cw + a cw + i/3 cw a ccvn (4) 

where a ccw and a cw represent the normalized energy am- 
plitudes of the CCW and CW modes, respectively; uj c 
assumes the same meaning as in Eq. (1), which de- 
notes the unperturbed resonance frequency of the WGM 
modes; Aw ccw and Aw cw are the resonance frequency 
shifts caused by the surface roughness to the CCW and 
CW modes, respectively; 7 CCW and 7 CW describe the corre- 
sponding scattering loss rates; /3 CCW is a parameter char- 
acterizing the coupling from the CCW to the CW modes 
and /3 CW is defined vice versa. From the fact that the 
CCW and CW modes only differ in their circulating di- 
rections, it is expected that 
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In fact, from Maxwell's equations, Aw ccw and /? C cw can 
be derived as [TT] 
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R _ uj c jAe(r)E* ccw (r)-E cw (r)d 3 r 
2 / e(r)\E ccw (r)\ 2 d 3 r 

where e(r) and Ae(r) correspond to the dielectric con- 
stant of a perfect microresonator and that of the surface 
roughness, respectively, and E ccw (r) (E CVJ (r)) is the elec- 
tric field of the CCW (CW) mode. Similar expressions for 
Alu cw and /3 CW exist with the exchange of the CCW and 
CW subscripts in Eqs. (6) and (7), and the relations in 
Eq. (5) become evident considering E ccvi (r) and E cw (r) 
are conjugate to each other (assuming e{r) is real). 

The scattering loss parameter 7 CCW (Tew) is obtained 
using the volume current method by computing the radi- 
ation power excited by the polarization current J ccw (r) = 
— iu] c Ae(r)E ccw (r). For example, the far-field electric 
field is given by [20] 



uj 2 c e lkar 



Ae(r')£ ccw (r') 



(8) 



(1 - rr)e- lknf - r ' d 3 r' 



where £rj and c are the permittivity and speed of light of 
free space, respectively; k is the wavenumber; and (r, 9, 
4>) are the spherical coordinates of the far-field position r, 
with (r, 9, </>) denoting the orthogonal unit vectors in the 
directions of increasing (r, 9, </>), respectively. In Eq. (8), 
we have adopted the exp (— iuj c t + ik$r) format for the 
outgoing light to be consistent with the convention used 
in Eqs. (3) and (4). The radiation power is calculated 
by integrating the Poynting vector over the sphere with 
radius r, and 7 CC w, according to its definition, is given by 
the power loss rate normalized by the mode energy as 



7ccw 



2je(r)\E ccw (r)\ 2 d 3 r 



(9) 



In addition, statistical information of the surface rough- 
ness can be inserted into Eq. (9), which leads to an en- 
semble average for 7 CCW (7 cw )[ll]. 

The eigenmodes of the coupled system are then ob- 
tained by solving the eigenvalue problem of Eqs. (3) and 
(4), which yields 
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where a± denote the two eigenmodes, whose resonance 
frequencies and scattering loss rates are given by 



w± = cj c +Au; ccw ± |/? ccw | 
7± = 7ccw 



(11) 
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FIG. 2. (a) Transmission response of a 20-/im-radius silicon 
microdisk which is pulley coupled to a 680-nm-wide access 
waveguide. Only the fundamental TE mode is phase matched 
and excited, (b)-(e) Zoom-in figures showing detailed line- 
shapes of the marked resonances. The resonant modes are 
over coupled so a higher extinction on resonance indicates a 
broader linewidth. 



Equation (12) predicts that a± should exhibit identical 
linewidths, while in experiments asymmetric lineshapes 
are often observed. For practical applications, j± are 
generally assumed to be different to fit the model to the 
experimental data |lllll8j. However, such an assumption 
contradicts with Eq. (12). The deficiency of this model 
can be further illustrated by an intriguing experimen- 
tal result depicted in Fig. 2, which shows the transmis- 
sion measurement of a 20-/mi-radius microdisk resonator 
fabricated on a silicon-on-insulator (SOI) wafer with a 
220-nm-thick silicon device layer |21j . By engineering 
the access waveguide geometry using the pulley coupling 
scheme [22], only the TE-polarized (electric field parallel 
to the device layer) first-order radial mode is excited, as 
confirmed by Fig. 2(a). Among this radial mode family 
(i.e., resonant modes with the same radial order and dif- 
ferent azimuthal orders), four resonances are picked out 
with zoom-in figures shown in Figs. 2(b)-(e), where sig- 
nificantly different lineshapes are observed. For Fig. 2(b), 
cj + w u)—, i.e., the mode splitting is negligible; for Figs. 
2(c)-(e), doublets appear and we have 7 + « 7_, 7 + < 7_, 
and 7+ > 7_ for each case. It seems difficult to explain 
the simultaneous occurrence of these features using the 
intuitive physics model, especially when the correspond- 
ing resonances belong to the same radial mode family in 
one specific microresonator. 

Summarizing the above discussions, we conclude that 
the two existing approaches for mode splitting and scat- 



tering loss have only achieved partial success. The 
independent-scatterer approach, based on collective con- 
tributions from each individual scatterer, works well 
when the number of scatterers is small and the scatterers 
are well separated with each other. On the other hand, 
the intuitive physics approach, based on a phenomeno- 
logical model for the mode splitting and the volume cur- 
rent method for the scattering loss, fits well for many- 
scatterer cases such as the surface roughness problem. 
However, both approaches have difficulties in providing 
correct scattering loss rates for the eigenmodes of the 
coupled system. In the independent-scatterer approach, 
lineshapes are predicted to be symmetric or asymmet- 
ric, but 7+ is always no more than 7_. In the intuitive 
physics approach, the model cannot predict asymmetric 
lineshapes (i.e., 7-1- 7^ 7_) in a self-consistent manner. 



III. A UNIFIED MODEL 

In this section, we will develop a model that is appli- 
cable to an arbitrary number of scatterers attached to 
the surface of a high-Q microresonator. The approach 
considers interactions among the CW and CCW modes 
and the free space continuum, with the coupling pro- 
vided by scatterers. The derivation here is similar to 
the independent-scatterer approach as in Refs. [5J [Tlj . 
but with a key modification that leads to distinct results. 
Moreover, we show conditions under which our model can 
be reduced to the results of the two existing approaches 
discussed in section II in their respective regimes. 

For clarity, we use the microdisk resonator as an ex- 
ample for the derivation, while the result is generally ap- 
plicable to any microresonator with a two-fold degener- 
acy in its resonance spectrum. For an isolated microdisk 
resonator (i.e., no external coupling), we can intuitively 
write down the following equation set: 



da r . 



dt 



dt 



= - \iu> c + 



Kq 



N 



m — cw,ccw 
N 



n—1 m— cw.ccw 



(13) 



(14) 



where a m and bj are the normalized energy amplitudes 
of the m (to = cw or ccw) WGM and the jth free space 
mode, respectively (lo c and u>j are their corresponding 
original resonance frequencies); kq is the intrinsic cavity 
loss without including the scattering loss; g n , m ,m' is a 
parameter describing the scattering of the m WGM to 
the same (to = to') or the counterpropagating (m ^ 
to') WGM mode induced by the nth scatterer; g n ,j,m 
is a similar parameter characterizing the nth-scatterer- 
induced scattering of the to WGM mode to the jth free 
space mode and g n ,m,j is defined vice versa. For now we 
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have used a discrete set of eigenmodes [bj] normalized in a 
finite but large enough volume to represent the free space 
continuum, and this restriction will be removed later. 

In this model, each scatterer is treated as a dipole. 
The electric field E excites the polarization of the nth 
scatterer as P — £oa n E n , where a n and E n are the po- 
larizability and the electric field at the location of the 
nth scatterer, respectively. The interaction of the polar- 
ization P with the electric field E is given by — P ■ E* 
[23] , with both the electric fields in P and E normalized 
to their corresponding mode energies. For example, the 
jth free space mode is expressed as 



1 



Ah a 



V^oVc 



n 



3 ' 



(15) 



where V c is the normalization volume of the free space 
modes; kj is the wave vector of the jth mode; and fij 
is the unit polarization vector. Similarly, the energy- 
normalized electric field of the m WGM has the following 
form: 



E m (r) = 



f(r) 



/[e(r)|/(r)| 2 d?r 



n r< 



(16) 



where we have explicitly written out the phase term 
exp (ik m x) with k m being the wavenumber of the m 
WGM along the mode circulating direction and x being 
the projection of r on that; f(r) accounts for the ampli- 
tude as well as the phase variation other than exp (ik m x); 
s(r) is the dielectric constant; and h m is the unit vector 
describing the polarization of the m WGM. To simplify 
the above expression, we can define a parameter V m as 



V m 



fe(r)\f(r)\ 2 d 3 r 



so 



and E m (r) can be alternatively expressed as 



E m {r) 



f(r) 



.e lk ™ x h rt 



(17) 



(18) 



Note that V m defined in Eq. (17) generally does not bear 
the unit of volume. However, we notice that f(r) is scal- 
able in Eq. (16). If we normalize the electric field of the 
WGM mode to that of a reference point, for instance, 
where the amplitude of the electric field is the maximum, 
f(r) can be interpreted as the relative field strength and 
V m defined above has the unit of volume. 

We now proceed to calculate the coupling coefficients 
9n,m,m', 9n, m , 3 , and g n ,j,m based on their definitions in 
Eqs. (13) and (14). Starting with Maxwell's equations, 
we have [20] 

, NN , x d 2 E(r,t) d 2 P(r,t) , . 

Vx(Vx£(r,f))+ M e(r) ^ 2 > = -» ^ \ (19) 



dt 2 



where u is the permeability of free space. For the m 
WGM mode, 

E(r,t) = a m (t)E m (r) = e-^'ta^tje 1 "^,^), 

(20) 



where we have separated the fast oscillating 
term exp (— iuj c t) with the slowly varying term 
a m (t) exp (iuj c t). Treating P(r,t) as a first-order per- 
turbation, Eq. (19) can be approximated as [2171 124j (also 
see discussions at the end of Appendix A) 

2±(a m (t)e iu ° t ) « iw c e*--*y P(r,t) ■ E* m (r)d 3 r. (21) 



P{r, t) consists of contributions from each scatterer as 



A' 



»(r,t) = eo5^a w (5^o m (t)B m (r) + 536 i (t)£? i (r)) 

3 

(22) 



n— 1 m 

x S(r - r„), 



where r n stands for the position of the nth scat- 
terer. Substituting the detailed expression of P(r, t) into 
Eq. (21) and comparing it to Eq. (13), we arrive at 



_ a n Uc\f(r n )\ A{k m ,-k m )x n 



2K 

_a n uj c f* 

9m ~ 2^y c 



(23) 



^ n - M & ■ ™M) , (24) 



where x n is the projection of the nth-scatterer's position 
r n along the WGM circulating direction. In Eq. (24), the 
dependence of the polarization of the m WGM mode n m 
on the position of the nth scatterer r n has been explicitly 
expressed, given that fi m is not necessarily a constant 
vector (for example, the dominant electric field for the 
TE-polarized WGM is E^, which is in the direction of 
4>). In deriving g n , m ,j, we have used the fact that only 
those free space modes that can resonate with the WGM 
modes (i.e., uj ~ oj c ) need to be considered, which will 
be proven soon. Taking a similar procedure for Eq. (14) 
leads us to 



(fij-h m (r n )J. (25) 



Note that g n , m ,j an d g n ,j.m obtained here are different 
from those obtained in Refs. [HI [33], where the derivation 
is based on the interaction among quantized fields and 
the exp(ikj-r n ) (exp(— ikyr n )) factor is missing in g n ,m,j 

(ffnj'.m)- 

With the knowledge of the coupling coefficients 
g n ,m,m< , 9n,m,j, an d g n ,j,m, w e are ready to solve Eqs. (13) 
and (14). Instead of studying the fast oscillating terms 
a m (t) and bj(t), it is more convenient to work with the 
slowly changing variables d m (t) = a m (t) exp (iui c t) and 
bj(t) = bj(t)exp(iuijt). We solve b 3 (t) in Eq. (14) in 
terms of a m {t) and insert it back into Eq. (13). After 
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some arrangements, we obtain 

N 

— = — — n_ it;) -4-7 



a m '(t) 

n—l ra' — cw,ccw 

/ 1 9n,m,j9n',j,m' 
ra' — cw,ccw n— 1 n' — 1 J 
i 

x / e^-^W-^a^t')^'. 



(26) 

Because the free space modes have a high mode density, 
we can replace the summation over the mode number j 
by an integral over the wave vector space as [2U] 



E 



Vc 



(27T) 



3 



TV TV OO 

J d <f> J sinOde J k 2 dk, (27) 



where (fc, t9, <j>) are the spherical coordinates of the wave 
vector fc, with (fc, 9, <p) denoting the orthogonal unit vec- 
tors in the directions of increasing (fc, (5, 0), respectively. 
rik describes two possible polarizations corresponding to 
fc, and hence can be in the direction of 8 or cj). 

Inserting the expressions of g n ,m,j arL d g n ,j,m into the 
last term on the right side of Eq. (26) , we arrive at 
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E EE 
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a n a n >uj 2 c f*(r n )f(r n ,) 
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4(27T) 3 V m C 3 
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sine 



/ w 2 e W»-»-»-n') 4,; / e l(w -^ )(t '- 4) (w(f/)di', 

— oo 

(28) 

where we have used the relation uj — kc to replace fc in 
the integral in Eq. (27). In the integration over t' , be- 
cause the microresonator has a high Q, it is reasonable to 
assume a m (t') varies sufficiently slowly over a few optical 
cycles so that it can be evaluated at the t = tf [5]. In 
consequence, 



j {u - UcW -t)- mi ^ dt , w ^ _ Wc) - m/ (t)) (2g) 



which indicates only the free space modes with resonance 
frequencies around uj c need to be considered. In Eq. (28), 
another simplification can be carried out for the summa- 
tion over the polarizations of the free space modes with 
the help of the following vector identity: 



(h-a){h-b) = (1- fcfc) -a (1 - fcfc) • b , (30) 



which can be easily verified for arbitrary vectors a and 
b. 

As a result, we have 



da m {t) 
dt 

with 

G 



= -ya m (<)+ ^ (iG m:r 



m'— cw.cew 



-)a m '(t), 
(31) 



= QL n u) c \f{r n )\ j( km ,-k m )x n 

m.m' — v t 7 & 



2V„ 



N N 



(32) 

n an>utf*( r n)f{rn>) j(-k m x n +k m ,x n ,) 



EE 

n—l n' — l 

[(1 - fcfc) • n m (r n ) ■ (1 - fcfc) • h m >(r n >) 



x e 



ikQk-(r n —r n / ) 



sin 8 d9dcj), 



(33) 



where fco = uj c /c is the wavenumber of light with an- 
gular frequency uj c in free space. The integral in T mtm i 
involves an integration over the spherical surface in the 
wave vector space, and its value depends on the polariza- 
tion of the WGM modes as well as the relative positions 
of scatterers. Hence, it is a geometric factor. 



A. Comparison with the independent-scatterer 
approach 

The geometric integral in T m m i given by Eq. (33) can 
be computed for any given r„ and r n >. For example, if 
the WGM mode is TM-polarized (magnetic field parallel 
to the device layer), n m = z (see Fig. 3(a)). In addition, 
we can choose the x axis to be in the direction of t*„ — tv . 
The geometric integral in Eq. (33) can then be simplified 
as 



I p&^O \ r n — r n ' I Sln & COS ( 



TV 

/87T 
sin 3 0J o {k o d n ,n' sin 9) d6 = —p{k d n ,n'), 



(34) 



where d n ^ n i = \r n — r n >\; Jo(x) is the Bessel function of 
the first kind of order zero; and p(x) is defined as 



p(x) 



7T 

J sin 3 6* Jo (x sin ( 



\d0. 



(35) 



n=e,, 



In deriving Eq. (34) , integral representations of the Bessel 
functions are used and mathematical details are left to 
Appendix A. From Eq. (34), we find that the geomet- 
ric integral in T m ^ m i is only a function of the separation 
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FIG. 3. (a) Illustration of the adopted coordinate system 
for the calculation of the geometric integral in Eq. (34): the z 
axis is perpendicular to the microresonator, and the x axis is 
chosen to be along the relative position of the two scatterers 
under consideration, (b) The solid line is the numerical result 
of p(x) defined in Eq. (35), and the dotted line is 3/2x, which 
corresponds to the envelop of p(x) when x is large. The exact 
expression of p(x) is given by Eq. (A7) in Appendix A. 



distance between scatterers; consequently, 



r , 

m.m,' 



N N 

EE 



a n an>u?:f*{r n )f(r, n 
6irV m c 3 



n=\n' = \ (36) 

In Fig. 3(b), numerical values of p(x) are evaluated. For 
x = 0, p(0) = 1; when the argument x is large, the 
envelope of p{x) decreases at a rate of which can 
be proved from a rigorous calculation (see Appendix A). 
In the independent-scatterer approach discussed in sec- 
tion II, only the n = n' terms (i.e., d n<n t — and 
p(0) = 1) are considered in the double summation over 
n and n' for T m ^ m i, and the n =/= n! terms are neglected, 
with the hope that the contribution from these terms 
is small if scatterers are well separated with each other 
|14j . From Fig. 3(b), we estimate that a reasonably large 
separation should be on the order of d n . n >/^o > 2.4 
{\p(kod nin >)\ < 0.1), with Ao corresponding to the free 
space wavelength (Ao = 2ttc/uj c ). Otherwise, the omis- 
sion of n 5^ n' terms can introduce significant errors and 
even lead to erroneous conclusions (such as 7+ < 7— , see 
discussions in section II). In section IV, we will examine 
one numerical example for two scatterers attached to the 
surface of a microresonator, where we show it is essen- 
tial to include p(x) for a complete understanding of the 
simulation results. 



B. Comparison with the intuitive physics approach 

When the number of scatterers is large, the forms of 
r TO)m ' given by Eq. (33) (or Eq. (36)) is not that conve- 
nient to work with. From Eq. (33), we notice that if we 



define 

S m ( 



JY 



47TvV m C 3 



E Q «/( r «)e ifcmXn nm(rn) 

n=l 

• (1 - kk)e' lkak - r - , 



^m,m' can be rewritten as 



S* (6,<j>) •S m /(6»,0) sin dd6d(f>. (38) 



The expression of S rn (9, cf>) in Eq. (37) is invariant if we 
scale f(r); hence, we can remove the restriction of E m (r) 
defined in Eq. (16) (or Eq. (18)) which requires it to be 
energy normalized, and extend the definition of S m {9, 4>) 
to arbitrary E m (r) as 



N 



4tt V U m c 3 ^ 

n—1 



^n-^m (^n) 



(39) 



(1 - kk)e 



— ikok-r n 



with U rn defined as 



U m = J £ (r)|£; m (r)| 2 d 3 r, 



(40) 



which corresponds to the energy of the m WGM mode. 
Also, any distribution of scatterers on the surface of a 
microresonator can be treated as a special case of surface 
roughness, with the dielectric perturbation given by 



Ae(r) 



N 



e a n S(r - r n ). 



(41) 



Combining Eqs. (39) and (41), a general form for 



is found as 



£0 



J Ae(r)E m (r) 



47r V U m c 3 J (42) 
■(l-kk)e- tkokr d 3 r. 

A direct comparison shows that, except a constant, 
S m (0, 4>) given by Eq. (42) in our model is equal to rE l c f w 
with £^' w given by Eq. (8) in the volume current method, 
if we identify (fc, 6, </>) (which are the coordinates of the 
wave vector fc) in Eq. (42) with (f, 9, 0) (which are the 
coordinates of the far- field position r) in Eq. (8). Fur- 
thermore, according to Eq. (31), T mtTn describes the scat- 
tering loss rate of the CCW (CW) WGM mode, similar 
to 7ccw in the phenomenological model in Eq. (3). Their 
difference is that T m m , given by Eq. (38), is a spheri- 
cal integration of \S m (9, cf))\ 2 in the wave vector space, 
and 7ccw, given by Eq. (9), is a spherical integration of 
|r£^| 2 in the coordinate space. From the mathemati- 
cal point of view, however, there is no difference in their 
detailed expressions and one can verify 7 CCW = T m>m . 



The analogy between our approach and the intuitive 
physics approach discussed in section II can be carried 
on further. With the help of Eq. (41), G m ^ m i defined 
in Eq. (32) can also be extended to the general case of 
surface roughness as 



components in the eigenmodes are not of equal weight 
(i.e., \r)\ is not necessarily equal to 1), in contrast to the 
result from Eq. (10), where the CW and CCW modes are 
equally weighted. Furthermore, from Eqs. (47)-(49), we 
have 



G ,, 



0J C 



As(r)E* m (r)-E m ,(r)d 3 r. (43) 



From here on, we will use a slightly different notation for 
the ease of comparison, that is, we use the subscripts m 
(> 0) and -m (< 0) to stand for the CCW and the CW 
WGM modes with the azimuthal order m, respectively. 
Comparing Eq. (43) to Eqs. (6) and (7), it is easy to rec- 
ognize Acj ccw = -G m>7n and /3 CCW = G m ^ m . Adding the 
already known relation 7 CCW = r m . m , the key difference 
between our model given by Eq. (31) and the phenomcno- 
logical model given by Eqs. (3) and (4) is that we have an 
additional coupling coefficient r m _ m in Eq. (31). After 
converting a m (t) back to a m (t) and setting kq to be zero 
(meaning only scattering loss is considered), Eq. (31) can 
be explicitly expressed as 



da. n 



dt 



da- 



ILU C — l(j m ,m i Z I a r, 



^G rn ^ m n J &_ m , (44) 



dt 



G, 1 m.m 
m,m H ^— I d- 



where we have used the facts G_ m _ r 



G, 



G- m>m = G*_ m , and r_ m , ro = T* m _ m , all of which are 



evident from their definitions (assuming Ae(r) is real). 
We also use the relation T- m - m = r m , m , which is not 
that obvious from Eq. (38), since from the expression of 
S m (M) gi ven b Y E q- ( 42 ), S m {9,dp) and S- m {9,4>) are 
not conjugate to each other. However, at current stage, 
we assume r_ m ._ m = Y miTn simply from the fact that 
the CW and CCW modes have the same scattering loss 
rate, and we will prove that in the next subsection. 
Solving the eigenmodes of Eqs. (44) and (45) yields 



1 



a.± 



C 



G m , m ± Re (t?(G* 



7± = r m>m T 21m (rj(G^_ m + -T* m ^ m )), 



(46) 

)), (47) 
(48) 



where 



(w± - lo c + G m , m )(7i 



-Re (G 



m,—m^ m,—m)i 



(50) 



(uj± — uj c + G m?m ) - (7± — r m)T7l ) /4 = 

Gfn—r 



J 2 /4. (51) 



For ui+, the first multiplying factor on the left side 
of Eq. (50) is no less than zero, and we con- 
clude that 7+ < 7_ if Re {G* m _ m T m! - m ) > and 
7+ > 7- if Rc (G* m _ m r m - m ) < 0. The case that 
Re (G* m _ m F m - m ) — is interesting, which corresponds 
to two possibilities, the first being the two eigenmodes 
have the same linewidth but different resonance fre- 
quencies (i.e., w+ > w_, 7+ = 7- = Y m ,m) and the 
second being the two eigenmodes have the same res- 
onance frequency but different linewidths (i.e., lj + = 
uj- = uj c — G m>m , 7+ 7^ 7_), depending on whether 



\G r . 



> 



,|/2 or \Gr, 



< 



t |/2 (see 



Eq. (51)). One trivial condition for Re (G 



is G„ 



, but since r r , 



r 

m, — m m,- 







generally is nonzero, 



the two resonances would have different scattering loss 
rates and can manifest themselves under different excita- 
tion conditions (i.e., weak and strong couplings). Thus, 
the two eigenmodes are not degenerate in the strict sense. 
However, both the experiment in Ref. [10] and the numer- 
ical study in Ref. [M] fail to observe the distinction of the 
scattering loss rates between the two eigenmodes when 
they overlap in the resonance frequency, mostly because 
only one transmission result with one particular excita- 
tion is available. In section IV, we will show numerical 
examples that clearly demonstrate when uj + = cj_, 7 + 
generally is not equal to 7_ . 

Next, we would like to see how the results given by 
Eqs. (46)-(48) can be reduced to those of the phenomeno- 
logical model derived in Eqs. (10)-(12). For dielectric 
perturbations in high-Q microresonators, usually (but 
not always) |G mj _ m | 3> |r m) _ m |/2. Under this con- 
dition, r\ given by Eq. (49) can be approximated as 



G 



i/G* 



, and Eqs. (46)-(48) are simplified as 



1 / Cm,, —771 

^ accwT W^\ acwl ' 



lj± m u> c — G„ 



± |G 



m,—m 1 1 



7± 



Re (G 



m,—m^ m,—mj 



G r , 



(52) 
(53) 
(54) 



V 



±4 



Gm. — 7 



IF 

2" 1 - m, — m 



^ m , — m ~^ 2 trt, — tn 



(49) 



The sign of 77 is chosen to ensure that uj + defined in 
Eq. (47) is no less than From Eq. (49), one im- 

mediate observation is that in general the CW and CCW 



Comparing Eqs. (52)-(54) to Eqs. (10)-(12), with the 
equalities of Aw ccw = -G„ vm and /3 CCW = G m _ m (see 
discussions following Eq. (43)), we find only 7± are dif- 
ferent between the two approaches. In fact, there is a 
simple explanation for the result of Eq. (54). The elec- 
tric fields corresponding to the eigenmodes a± are the 
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eigenvectors of Eqs. (44) and (45), which can be solved 
as 



E ± (r) 



1 t G„ 



v2 ^ |G m ,_ 



(55) 



The associated scattering loss rate can then be calcu- 
lated based on the volume current method. The far-field 
electric fields corresponding to E± (r) have the same lin- 
ear combinations as in Eq. (55) by those of the CW and 
CCW modes. Thus, 



E f ±(r)\ 



E, 



far 



. .a Re(G 
(r)\ =F- 



m. — m rn 



E^*(r)E^ m (r)) 



IG, 



(56) 

The scattering loss rate involves an integration of 
E± 1 (r) over the sphere with radius of r as in Eq. (9) . 
Using the equivalence we have established between 
S m (8, <j>) and r£^ r (see discussions following Eq. (42)) as 
well as the expression of T mtm i given by Eq. (38), we find 

if 1 2 

the spherical integration of LE^ r (r) is equal to T my7n 
and the spherical integration of E^* (r)E^ m {r) is equal 
to r m> _ m , and Eq. (54) becomes apparent. 



C. Formulation in the Fourier domain 

In this section, we will derive a useful formulation for 
the calculation of parameters needed to obtain mode 
splitting and scattering loss, i.e., G m , m , G m _ m , Y m ^ m 
, and T m ._ m in Eqs. (44) and (45). We still use the 
microdisk resonator as an example, and we further as- 
sume that the scatterers are uniform along the microdisk 
slab thickness (this contains the Rayleigh-scatterer case, 
for which the detailed shape of scatterers is not impor- 
tant). As a result, the surface roughness Ae(r) defined 
in Eq. (41) can be assumed to have the following form: 



Ae{r) = £ 5(r - R)Ae r (<f>)rect(z/h), 



(57) 



where R is the radius of the disk; (f> is the azimuth, as 
shown in Fig. 3(a); Ae r ((f>) characterizes the dielectric 
perturbations along the periphery of the microdisk res- 
onator; h is the slab thickness; and rect(x) stands for 
the rectangular function [35]. Because of the inherent 
periodic boundary condition, As r ((f>) can be expanded in 
terms of periodic harmonics along the microdisk periph- 
ery as [IB] 



(58) 



where F(k n ) is the Fourier component of Ae r (<f>) with 
k n = n/R (n = 0, ±1, ±2, . . . ). The Fourier transform 
of Eq. (58) gives 



2tt 



F(k n ) = J Ae r {<t>)e- in * d<j>. 
o 



(59) 



For dielectric perturbations, Ae r (4>) is real, which yields 
F(k n ) = F*(k- n ). Inserting the form of Ae r ((f>) given by 
Eq. (58) into Eq. (43), G m m and G m ._ m are obtained as 



with 



G m ,m — goF(ko), 

III } 



Enu c Rh , - ,2 
9o = \E m (R,0)\ , 



(60) 
(61) 



(62) 



where \E m (R, 0)| is the amplitude of the electric field 
E m {r) at the surface (r = R, <p = 0) after averaging 
along the slab thickness. 

To obtain T m ^ m and T m — m , we first calculate S m (8, (j)) 
based on Eq. (42). For the TM-polarized WGM modes, 
h m = z , and Eq. (42) can be simplified as [16] 



x / e 



a i(m— n)<f>' — ikoH sin 9 cos (0 - 



= v / G "(- sin 90) F*(k m+n )(ie^r n J n (k Q Rsme), 

n 

(63) 

where the integral representation of J n (x) is used (see 
Appendix A) and Go is defined as 

Likewise, we have 

S^ m {6,<j>) = ^Co(-sm80)Y,F(k m+n )(ie- l *y n 

n 

x J n (k Rsin9). (65) 

Substituting S m (9, 4>) and S_ m (9, <j>) into Eq. (38), r m m 
and r m ,_ m are found as 

= 2 7 rG ^|F(A: m+ , i )| 2 

n 

x J Jl (k a R sin 6) sin 3 6d6, (66) 
o 

and 

L m,—m ^ 

7TG ^F(fc )F(k m - n ) 
n 

7T 

x J Jl ( k aR sin 6) sin 3 9d6. (67) 
o 

One can verify that if we replace m by — m in Eq. (66), 
after some arrangements (replacing the summation in- 
dex n by — n and using _F(fc_„) = F*(k n ), J_„(x) = 
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(— 1)™ J n (x)), it will lead to the same expression. This 
proves r m m = r_ mi _ OT , which we have already used in 
Eqs. (44) and (45). ' 

For the TE-polarized WGM modes, the dominant elec- 
tric field is E^, which is in the direction of cj>. Equation 
(42) then becomes [TB] 



S m (9,<p)=y/Co4>^-y2F*(k n ) 

Z7T * — ' 

n 

7T 

x J e i(m-n)(f>'-ikoR Bin 6 cos ($'-<$) cos (^/ _ 

—TV 

n 

x ^ J n -x (koR sin 9) — J n +i (koR sin 8)j . 
Similarly, we have 



(68) 



-m\V,ip) =- 



(69) 



x (j n -i(koRsia9) — J n+ i(k Rsi 
It follows from Eq. (38) that 



7T 

: J (j n -i(kaRsin6) - J n+ i(k Rsm 



smt 



r - nCa 



sin 9 d9, 
(70) 



and 



m. — m 



(71) 



IV. APPLICATIONS 



To demonstrate the applicability of the developed 
model, in this section we will apply it to three different 
examples. The first two examples deal with one and two 
scatterers, respectively, and in the third example we con- 
sider the fabrication-induced surface roughness present in 
high-Q microdisk resonators, which corresponds to thou- 
sands of small scatterers. Numerical and experimental 
evidences are provided to support the derived theoretical 
results. 



A. Single scatterer 

For the single-scatterer case, Ae r {(j)) = aS(4> — (j)o), 
where a is a constant and 4>o is the position of the scat- 
terer. From Eq. (59), we have F(k n ) = aexp(— intfio). 
It follows from Eqs. (60) and (61) that G m . m — goa, 
and G m ._ m = goa exp(— i2m<fio). For the TM-polarized 
WGM mode, T m m is given by Eq. (66) as 



7T 

T m , m = 27rC a 2 J JlikoRsinO) sm 3 6d6, 



(72) 



where the following Bessel identity has been used (see 
Appendix A): 



Similarly, T r , 



J24(x) = l. (73) 

n 

is obtained from Eq. (67) as 



(74) 



r\ defined in Eq. (49) can then be calculated to be 
1] = exp (— z2m</>o), which indicates that the CW and 
CCW modes are equally weighted in the eigenmodes a±. 
Moreover, from Eqs. (47) and (48), we have 



with 



uj± — u c — Go ± Go, 
1± = r T r , 



2 

G = goa; T ee -7rG a . 



(75) 
(76) 



(77) 



Equations (75) and (76) reproduce the familiar results 
for the single-scatterer example, which are independent 
of the position of the scatterer (i.e., 4>o)- This stems 
from the fact that physically measurable scalar variables 
should be invariant with respect to the choice of the co- 
ordinate origin. Therefore, one can simply take <po = 
and arrive at the same results in Eqs. (75) and (76). 

The example mentioned in Fig. 1(b), which is equiv- 
alent to a negative-dielectric-constant scatterer on the 
surface of a larger-radius microresonator, can also be an- 
alyzed. Assuming the scatterers have uniformly covered 
the surface of the microresonator except for a vacancy at 
6o = 0, we have Ae r (<j>) = ^{4> ~ 4>i) ~ &{4>) (f° r sim- 
plicity, we neglect a constant here, i.e., a = 1), where fa 
is a set of angles representing the locations of these small 
scatterers (with the vacancy filled too, since we subtract 
it in the second term of Ae r ((/>)), which uniformly fall in 
the range of (0, 27r). The first term in Ae r (<j)) only con- 
tributes to F(ko), and we have F(k n ) = —1 for « / 0. 
As a result, G m - m = — G ■ Though F(k ) is a large 
number, it does not contribute to the scattering process, 
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given that physically it corresponds to a uniform thin 
layer of dielectrics. Mathematically, from Eqs. (66) and 
(67), we find that the weight coefficient for F(k ) is pro- 
portional to an integral of J m (fco-Rsin#) . One impor- 
tant property of J n (x) is that its value is only significant 
when \n\ < \x\. Because m > k$R (m — koRn c g, with 
n e g being the effective index of the WGM mode), the 
contribution of F(ko) to r m m and r m ,_ m is negligible. 
Hence, T mtm and T m ^ m are the same as those of the 
single-scatterer case given by Eqs. (72) and (74), respec- 
tively. This leads to Re (G^ _ m T m _ m ) = -G T < 0, 
and 7+ > 7-, as expected. 

B. Two scatterers 

We start with two identical scatterers. Since we have 
argued that only the relative positions of scatterers are 
important and the choice of the azimuthal origin can be 
arbitrary, we can take Ae r ((f)) = 5{4> — 0o) + <$(</> + M 
(again, we have omitted a constant in this expression). 
As a result, F(k n ) = 2cosn</> , G TO , m = 2Go, and 
G m> - m = 2G Q cos2mcf) . For the TM-polarized WGM 
mode, it follows from Eq. (66) that 

T m ,m = 87tG ^ cos2 (( TO + n )^o) 

n 

x J J^(k Rsm9) sin 3 6 d9. (78) 
o 

The above result can be simplified using the following 
Bessel identities (see Appendix A): 

Jn (koR sin 6>) sin 2n<J) = 0, (79) 

n 

J^(koRsm6) cos 2n(f>a = Jo(2k Rsm9 sin^ ), (80) 

n 

and we get 

T m , m = 2r (1 + p(2k R sin 4> ) cos 2m<j> ) , (81) 

where p(x) defined in Eq. (35) is used. Similarly, from 
Eq. (67), 

r m _ m = 2r (p(2k R sin <p ) + cos 2m0 o ) ■ (82) 

Equations (81) and (82) can also be derived from 
Eq. (36), which is more convenient for this case. Since 
G m - m and r m- _ m are both real, from Eq. (49), r\ = ±1. 
By our convention, the choice of 77 is to ensure w+ is no 
less than w_. Thus, according to Eq. (47), 77 takes +1 
when G m _ m > and takes —1 when G m _ m < 0, and 
can take either 1 or —1 when G m ._ m = , at which point 
the two eigenmodes share the same resonance frequency. 
Equation (47) then becomes 

uj ± =Lu c -2G Q ±2G Q \cos2m(f>o\. (83) 



In addition, j± is given by Eq. (48) as 

7± *^m,m "F ?]^rn,—m 

—2T (1 + p(2k n R sin </> ) cos 2m</> ) 

=F 2r (| cos2m0 o | + sign(cos2TO(/) )p(2fc i? sin </>□)) > 

(84) 

where in Eq. (84) we have substituted r\ by the sign func- 
tion of cos2m</>o, which only differs with rj at the zeros 
of cos 2mcf>Q . There are two reasons that allow us to do 
that. First, at the zeros of cos2mcf>o, the definition of 
7± becomes ambiguous, since the subscripts ± are only 
used to distinguish the resonance frequencies of the two 
eigenmodes. Second, as we shall show below, the zeros 
of G m - m (oc cos2?n</>o) are singular points of 7±, where 
the behavior of 7± can only be studied by infinitely ap- 
proaching these points. 

As has been discussed in the comparison with the 
independent-scatterer approach, when the two scatterers 
are separated at large distances, p{x) can be neglected in 
Eq. (84) and we have 

7± «2r o T2r o |cos2m0 o |, (85) 

which is identical with Eq. (2) from the independent- 
scatterer model (N — 2). However, when the two scat- 
terers are close to each other, p{x) has to be considered 
in Eq. (84) . One special example is that the two scatter- 
ers overlap with each other (i.e., <fio = and p(0) = 1), 
which can be treated as a single-scatterer case. Equation 
(84) then predicts 7+ and 7_ to be and 8r , respec- 
tively. In contrast, Eq. (2) provides 7+ and 7_ to be 
and 4r , respectively. Using the result of Eq. (76) for 
the single scatterer (7+ = and 7_ = 2Fo), we find 
that Eq. (84) is accurate and Eq. (2) only predicts half 
of the exact number for 7_(as seen from Eq. (77), the 
scattering loss is proportional to the square of the di- 
electric perturbation so should be four times bigger if 
the dielectric perturbation doubles). Another interesting 
observation is that when approaching the zeros of the 
cos2m</>o, p(2koRsm(j)o) is generally nonzero, and 

7± w 2r T 2r o sign(cos2m0 o )p(2fco-Rsm(?!>o)- (86) 

If we sweep 4>o continuously, each time cos 2m<f>o crosses 
its zero points, its sign will change and there will be 
abrupt changes in 7 + and 7_, indicating the zeros of 
G m ._ m (cx cos2m</>o) are singular points of 7 + and 7_. 
Moreover, since the relation between 7+ and 7_ are re- 
versed when passing the zeros of cos2m</>oi it is always 
possible to observe 7+ > 7_ in the neighborhood of 
G m -m — (as long as p(x) is not negligible). 

To verify the derived theoretical results, we perform a 
numerical investigation for a two-scatterer example using 
an in-house two-dimensional (2-D) microresonator mode 
solver implemented in the COMSOL environment |26j . 
Details of the implementation are provided in Appendix 
B. The inset of Fig. (4) illustrates the studied structure, 
which consists of two 10-nm-radius scatterers attached to 
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FIG. 4. Simulation results of Wdiff and 7diff, which are de- 
fined by Eqs. (87) and (88), respectively, for two 10-nm-radius 
scatterers attached to the surface of a 2-/zm-radius microdisk 
resonator as illustrated by the inset. The refractive index 
of the microdisk is 2.829 (obtained using the effective index 
method for a 220-nm-thick silicon layer), and the refractive 
index of the scatterers is twice as big (i.e., n sca tt C rcr = 5.658) 
to make the scattering effect significant, m is the azimuthal 
order of the WGM mode (TM polarized), which is 19. uj c is 
obtained from the simulation for an ideal microdisk resonator 
without any scatterers as 1.2282723el5 rad/s, and Go and To 
are obtained from the single-scatterer simulation as 5.88el0 
rad/s and 9.3e8 rad/s, respectively. 



the surface of a 2-/im-radius microdisk resonator. We fix 
the position of one scatterer and sweep the position of the 
other. The complex eigenfrequencies of the coupled sys- 
tem are computed by the mode solver, offering both the 
resonance frequencies and the scattering loss rates for the 
two eigenmodes. In Fig. (4), two normalized parameters 
Wdiff and 7diff are plotted, which are defined as 



l^diff 



7diff : 



4G 

7+ - 7- 

4r 



(87) 



From Eqs. (83) and (84), our theoretical model predicts 



7diff 



^diff = | cosm0|, (89) 
| cosm<t>\ — sign(cosm0)p(0), (90) 



where the angular separation between the two scatterers 
are <fi = 2<^o as shown in the inset of Fig. (4) . Comparing 
Fig. (4) to Eqs. (89) and (90), we find Wdiff agrees with 
Eq. (89) well; and 7diff indeed changes sign when passing 
the zeros of Wdiff- We also notice that the magnitude 
of 7diff can be less than —1; especially, it approaches to 
—2 when the two scatterers are close to each other, as 
expected from the single-scatterer result. In Fig. (5), we 
plot two additional normalized parameters w sum and 7 SU m 



FIG. 5. Numerical results of w sum and 7 sum , which are de- 
fined by Eqs. (91) and (92), respectively. The red triangles 
corresponds to 7 sum directly from simulation, while the black 
line corresponds to 7s Um obtain by extracting p(4>) from the 
simulation result of y^iB first (using Eq. (90)) and then com- 
puting the numerical values of Eq. (94). 



defined as 



7s 



(91) 
(92) 



4G 

. 7+ +7- 

4r ' 

which are shown by the blue solid line and the red triangle 
marks, respectively. According to Eqs. (83) and (84), 



-1, 



-'sum 

7sum = l+p(</>) cosm^. 



(93) 
(94) 



As observed from Fig. (5), w sum has a few percent fluctu- 
ations around the theoretical value (i.e., —1), largely aris- 
ing from the limited positioning resolution of the moving 
scatterer when we sweep it along the perimeter of the 
microdisk (1 nm in the COMSOL environment). With 
the help of Eq. (90), we could extract p(<j>) from the nu- 
merical result of 7diff shown in Fig. (4), and the result 
is depicted by the red triangles in Fig. 6. Moreover, us- 
ing the obtained p(<p) , 7 sum could be computed based on 
Eq. (94). The result, which is shown by the black solid 
line in Fig. 5, agrees with the one from direct simulation 
(red triangle marks) well, implying that our theoretical 
model is self-consistent. One may notice that p(<fi) shown 
in Fig. 6 is different from the one plotted in Fig. 3(b). 
This is because p((f>) shown in Fig. 3(b) is for the three- 
dimensional (3-D) case, while our simulation considers a 
2-D model. The essential difference can be traced back 
to the difference in the free-space Green's function |20) . 
Employing the 2-D free-space Green's function and fol- 
lowing a similar procedure as in the 3-D case(see the end 
of Appendix A), we obtain 



p{4>) = Jo(kod) = Jo ( 2fcoi?sin 



(95) 
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FIG. 6. The red triangles corresponds to p(<fi) extracted from 
the simulation result of fdiff shown in Fig. 4 based on Eq. (90). 
The black solid line is the numerical result of Eq. (95) , which 
is the theoretical prediction of p(<f>) using the approximate 2- 
D Green's function. The blue dashed line is the numerical 
result of Eq. (Bll) (Appendix B), which is the theoretical 
prediction of p(4>) using the accurate 2-D Green's function. 



where d is the distance between the two scatterers. This 
result can also be expected from Eq. (35), by taking the 
inclination coordinate 9 — n/2 and skipping the inte- 
gration over 9. In Fig. 6, we have plotted the predicted 
p(<fi) given by Eq. (95) by the black solid line, which 
agrees with the one extracted from the numerical sim- 
ulation (shown by red triangles) reasonably well. The 
deviation there arises from two facts. First, we have cer- 
tain positioning error when sweeping the scatterer in the 
simulation, as already mentioned for w sum in Fig. 5. Sec- 
ond, an accurate p(4>) requires taking the effect of the 
microdisk resonator to the free-space Green's functions 
into account, which has been omitted in Eq. (95) (or 
Eq. (35)) (see the end of Appendix A for more discus- 
sions). In the Appendix B, a brief derivation is provided 
for the accurate calculation of p{4>) in the 2-D space, and 
the result is shown by the blue dashed line in Fig. 6, which 
agrees with the simulation result well. From the asymp- 
totic behavior of Jo(x) (Appendix A), one notice that in 
the 2-D case, the magnitude of p{x) decreases with the 
separation distance d of the two scatterers as 1/ ^kgd, in- 
stead of l/k^d as in the 3-D case (Fig. 3(b)). Therefore, 
a larger separation distance is required for 2-D models to 
neglect the effect of p((j>) (d/X > 8 for \p(4>)\ < 0.1). 

Finally, we would like to mention that if the two scat- 
terers are not identical, by properly choosing the origin 
of the azimuth, we can still make G m - m (which only 
depends on Fik^m)) real, but r mj _ m (which depends on 
multiple terms of F(k n ) as given by Eq. (67)) generally 
will be complex. From Eq. (49), 77 does not necessarily 
have a magnitude of 1, which means the CW and CCW 
modes are not equally weighted in the eigenmodes. To 
observe a significant deviation of \r]\ from 1, |r m _ m | has 
to be close to |G m ._ m |, which could only happen when 



Gm-m is near its zero points, since for dielectric scatter- 
ers Go 3> To (as an example, for the scatterers studied 
in Fig. 4, Go/r = 63). Thus, our model offers a simple 
explanation for the inequality of the CW and CCW com- 
ponents in the composition of the eigenmodes, which is 
studied in detail in Ref. 115 . 



C. Fabrication- induced surface roughness 

Here we will examine a different example, i.e., the side- 
wall roughness caused by the imperfect fabrication of 
microresonators, which corresponds to numerous small 
scatterers on the surface. The distribution of these scat- 
terers is random and usually follows a stationary statistic 
as [27] 

< Ar(x)A(x') >= a 2 exp ( J x ~ x l ) a ( 96 ) 



i\ > _ 2 / \ x x '\ 2 



< Ar(x)A(x') >= (7 exp ( 



Li 



), (97) 



where Ar(x) is the sidewall roughness along the wave 
propagation direction x; <> stands for the ensemble av- 
erage; a is the roughness standard deviation; and L c 
is the correlation length. For the microdisk resonator, 
the dielectric perturbation function Ae r (</>) is related to 
the sidewall roughness as Ae r (cf>) = 8n 2 Ar(R(f>), where 



Sn 2 = 



Uq , with rid and no being the refractive in- 



dices of the microresonator and the surrounding medium 
(air in our case), respectively. Substituting Eq. (58) into 
Eqs. (96) and (97), we obtain 

< F(k n )F* (k m ) >= t^Z a l L ^ n - m )> (98) 



and 

<F(k n )F*(k m ) >-- 



R{l + {k n L c f) 



2-7T2 (Sn 2 ) 2 a 2 L c k n L c 2 



R 

x 5(n — m), 



expH^pn 



(99) 



respectively. The Kronecker's delta function in both 
Eqs. (98) and (99) implies that [F(k n )} (n > 0) are 
statistically independent random variables (remember 
F(k n ) = F*(k- n )). Specifically, each resonator is one 
possible realization of [F(k n )], and Eqs. (98) and (99) are 
valid when the ensemble average is performed for many 
independently fabricated resonators under the same con- 
dition. 

The independence of \F(k n )] (n > 0) provides key 
insights to the understanding of the properties of mi- 
crodisk resonators such as the one shown in Fig. 2. Us- 
ing the results of Eqs. (52)- (54) under the assumption 



of G„ 



> |r„ 



J/2, we find the mode splitting 



for the azimuthal order m is proportional to |G mj _ r 
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and therefore |-F(fc2 m )| (Eqs. (53) and (61)). Because of 
the independence of [F(k n )], different azimuthal orders 
can have independent mode splittings. Strong variations 
of mode splitting over azimuthal orders for an individ- 
ual microresonator are thus possible if the correspond- 
ing [F(k n )] have strong amplitude variations with the in- 
dex n. Moreover, [F(k n )] should also have independent 
phase variations. For example, if [F(k n )] are all posi- 
tive numbers, then the correlation function in Eqs. (98) 
and (99) will also be positive, contradicting with the re- 
sults there. The variations of the phase and amplitude of 
[F(k n )] (with the index n) can lead to variations of scat- 
tering losses for different azimuthal orders. To see this, 
we use the TM polarization as an example. In Eqs. (66) 
and (67), T m m and T m _ m for the azimuthal order m 
is summed over terms of [F(k n )] with indices around 
m, weighted by coefficients as integrals of J n (koRs'm6). 
As already mentioned, J n (x) only has significant values 
when \n\ < Therefore, for azimuthal order m, the 
summation in Eqs. (66) and (67) for T mim and r m ._ m 
only contain limited terms of [F(k n )], with contributing 
indices in the range of (m — k R, m + k R) (see Fig. 7 for 
an illustration). This has two consequences. First, r m m 
and r mi _ m will show some dependencies on m because of 
this local summation cannot average out the variations 
among [F(k n )]. Second, whether 7 + > 7_ or 7 + < 7_ is 
determined by Re (G* m _ m r m ,_ m ), which is given by 

Rc (G^_ m r m ,_ m ) cx Re (F*(k 2m )T m ^ m ). (100) 

Since m > k^R, F(k2m) does not correlate with r m ._ m 
(see Fig. 7). If the phase of F(k2 m ) can vary randomly 
within (0, 2tt), then Re {G* m ^_ m T m ^ rn ) can be either pos- 
itive or negative. Furthermore, because of the indepen- 
dence [F(k n )], Re (G* m _ m T m - m ) will also be indepen- 
dent for different azimuthal orders. For that reason, 
one can observe simultaneous occurrence of different line- 
shapes among the same radial mode family in an individ- 
ual microresonator, as shown in Fig. 2. The above discus- 
sions also apply to the TE polarization, since the weight 
coefficients in Eqs. (70) and (71) have similar properties 
as those of Eqs. (66) and (67). 

The next question is what the amplitude and phase 
variations of [F(k n )\ with the index n for an individual 
microresonator. One method is to generate scatterer dis- 
tributions following the statistical rule given by Eq. (96) 
or (97), and calculate [F(k n )} according to Eq. (59). We 
have constructed one such scatterer distribution in Ap- 
pendix C, which satisfies Eq. (96). Numerical simulations 
there indicate that the amplitude of each F(k n ) follows 
a Rayleigh distribution and its phase follows a uniform 
distribution within (0, 2-k). A more convenient way is to 
assume some simple statistical models for \F(k n )]. For 
example, for the surface roughness that can be approxi- 
mated by Eq. (97) , in Ref. [TB] we have assumed a Gaus- 
sian distribution for the amplitude and uniform distribu- 




Azimuthal order 



FIG. 7. The grey thin bars depict [F(k n )], and the two red 
thick bars correspond to F{k m ) and F{kim) for the azimuthal 
order m under consideration. Ym,m and r mj _ m , given by 
Eqs. (66) and (67), respectively, are limited sums of [F(k n )] 
with indices around m, with the weight coefficients given by 
the black solid line. On the other hand, G m ,-m is propor- 
tional to F(k2m), and therefore has no overlap with r m , m or 
r 

tion for the phase of F(k n ) as 

x (cos a + sin a N n (0, 1)) exp (i2nU n (0, 1)), 

where [N n (0, 1)] are independent random variables with a 
normal distribution with a zero mean and a unit variance, 
and [U n (0, 1)] are independent random variables with a 
uniform distribution in the interval (0, 1). The parameter 
a is introduced to account for the amplitude variations 
of [F(k n )]. The independent uniformly distributed phase 
terms [exp (i2irU n (0, 1))] will ensure the independence of 
[F(fc„)], as required by Eq. (99). With the generated 
[-F'(^n)]: G mjm , G mj _ m , T m m , and T m _ m are obtained 
for each azimuthal order m, and the mode splitting and 
scattering loss rates are calculated from Eqs. (47)-(49). 

Figure 8 shows one simulation example for a 10-^tm- 
radius silicon microdisk resonator. We have defined 
dimensionless quality factors for the scattering loss as 
Q ss ,± = ^c/l± and for the mode splitting as Qp = 
uj c /(oj + — w_). We also define a scattering Q for the 
CW (CCW) mode as Q SSt t = ^ c /r m , m , though it is not 
directly measurable since the CW and CCW modes are 
no longer the eigenmodes of the system. However, from 
Eq. (48), we have 

2Qj, Q.:. . ■ G • (102) 

In Fig. 8(a), we have plotted Q ss ,± for the two eigen- 
modes in the blue dotted line with circles and red dotted 
line with squares, respectively; Q ss ,t is plotted in the 
black dashed line while Qp is shown by the green dash- 
dotted line with crosses. We have chosen the parameters 
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FIG. 8. (a) Simulated mode splitting and scattering loss for 
a 10-/im-radius silicon microdisk resonator with \F(k n )] gen- 
erated from Eq. (101). The parameters used are a = 2 nm, 
L c — 160 nm, and a — 0.5n. \E m (R, 0) is evaluated by 
obtaining the fundamental TE WGM mode from 3-D finite 
element method (FEM) simulations 18 , followed by a subse- 
quent averaging along the vertical dimension of the microdisk 
slab. The dependence of \E m (R, 0)| on m (or the wavelength) 
is small and can be neglected, (b)-(e) Transmission responses 
of the marked resonant modes with a coupling Q of 600,000 
(over-coupled). Therefore, similar to Fig. 2, a higher extinc- 
tion on resonance indicates a broader linewidth. 



in Eq. (101) to generate close results to the experimen- 
tal ones as shown in Fig. 2, with a = 2 nm, L c = 160 
nm, and a = 0.5ir. As observed from Fig. 8(a), indeed 
we have strong variations of mode splitting for different 
azimuthal orders. Furthermore, as shown by the four 
zoom-in figures in Figs. 8(b)-(e), different lineshapes can 
be observed for different azimuthal orders in the trans- 
mission measurement, similar to the experimental results 
shown in Figs. 2(b)-(e). In particular, in Fig. 8(e), only 
a single resonance is observed because the mode split- 
ting is negligible. However, the resonance actually corre- 
sponds to two eigenmodes with different scattering loss 
rates, and will show different responses under different 
coupling conditions. Another point worth mentioning is 
that the scattering Qs of the both eigenmodes exhibit 
more than 30% variations over azimuthal orders, and in 
a scattering-loss limited microdisk resonator, such varia- 
tions will be transferred to the intrinsic Qs [T6HT8] . Fi- 
nally, because of the random nature of Eq. (101), Fig. 8 
is just one possible result, and different simulation runs 
will generate close but not exactly the same outcomes. 
This in fact closely mimics the real fabrication, which 
produces resonators with comparable but rarely identi- 
cal performances. 



V. CONCLUSIONS 

In summary, we have developed a unified model that 
applies to an arbitrary number of scatterers, which 
provides a comprehensive understanding on the mode 
splitting and scattering loss in high-Q WGM microres- 
onators. Compared with the independent-scatterer ap- 
proach which is commonly used for the a-few-scatterer 
scenario, our work reveals that the independent-scatterer 
model has neglected the interference terms from differ- 
ent scatterers, whose effect decreases with the separa- 
tion distance d as 1/kod for 3-D cases and as l/y/k^d for 
2-D cases. Thus, the independent-scatterer model only 
works when scatterers are well separated. Compared 
with the intuitive physics approach which is developed 
for the many-scatterer scenario, we have derived an ad- 
ditional coupling term between the CW and CCW modes 
(i.e., r mj _ m ) that has been missing in the phenomenolog- 
ical model used by the intuitive physics approach. This 
modification leads to the prediction of asymmetric line- 
shapes in a self-consistent manner. Moreover, combined 
with numerical studies and experimental results, the uni- 
fied model has provided many new understandings on the 
mode splitting and scattering loss in high-Q WGM mi- 
croresonators. For example, we prove that the intuitive 
belief that 7+ < 7_ is not generally true, and counter ex- 
amples can even be found for two scatterers attached to 
the surface of WGM microresonators. Our work also un- 
veils that when mode splitting disappears, the scattering 
loss rates of the two eigenmodes are generally different, 
and 7+ and 7_ become singular at these points. For the 
fabrication-induced surface roughness present in high-Q 
microresonators, the stationary distribution of numerous 
small scatterers results in independent mode splitting for 
different azimuthal orders. The scattering loss rate of 
each eigenmode will also show more than 30% variations 
among the radial mode family, which explains the ob- 
served variations of intrinsic Qs in scattering-loss-limited 
microresonators [TBI E! • We believe such a unified ap- 
proach does not only fill the gap for the existing theo- 
retical works on the mode splitting and scattering loss 
in high-Q WGM microresonators, but also will play an 
indispensable role for the practical applications of these 
phenomenons to produce the most accurate results. 
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Appendix A: Mathematical Formulas 

In this part, we list the mathematical formulas that 
are used in the paper, mainly for the Bessel functions. 
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Brief derivations are provided for some equations which 
are not easily found in mathematical handbooks. 

J n (x), the Bessel function of the first kind, has the 
following integral representation [25] : 



Jn{x) = / e^ xsinT - nT Ur. (Al) 

Z7T 



Substituting r = tt/2 — r', we obtain another represen- 
tation of J n (x) as 

7T 

J n ( X ) = ^-i~ n I e^osrinr)^ ( A2 ) 

2n I 



we have 



,i2n<j> a t2 



J 2 (x) 



2tt 



i(x sin T\ -\-x sin T2 



>S(ti + r 2 - 2(/> ) dr\dri 



(A10) 



Changing the integration variables from n and T2 to (ti+ 
t 2 )/2 and (ri — r 2 )/2 would lead us to 

^y 2 "*°j2(a;) = J (2a; sin O ). (All) 

n 

In particular, if we let (j>o = 0> Eq. (All) becomes 

E J n^) = L ( Al2 ) 



When x 3> n 2 , we have the following asymptotic forms 



J n (x) 
H^{x) 



2 nf 1 

cos (a; - — --), A3) 

7TX 2 4 

2 mr 7r , 

exp(t(a?- — -- , A4 

7TCC 2 4 



where Hn (x) is the Hankel function of the first kind. 

To obtain p(a:) defined in Eq. (35), we employ Sonine's 
first integral as [2Ej 



tt/2 

/ 



Jo (2 sin0) sin6>cos^ +i 6>d<9 



2 v T(v + l)J u+1 {x) 



r y+l 



where Y(v + 1) is the gamma function |25j . p(x) can be 
evaluated by taking v in Eq. (A5) to be — | and | and 
subtracting these two terms, which results in 



(A6) 



Equation (A6) can be further simplified with the help of 
the spherical Bessel functions [25] . and we obtain 



. . 3 / sin x 



cos x sin x 



(A7) 



When x is large, the leading term of p(x) is the sine 
function and its envelop varies as 3/2x. 

Equations (73), (79), and (80) can be proved by consid- 
ering the following series, which can be calculated using 
the integral representation of J n (x) given by Eq. (Al) as 



4tt 2 



^ ' e i(xsinT 1 +xsinT2-nT 1 -nT 2 +2n<fio) ^.^^ 



(A8) 



Using the following identity: 



(A9) 



Regarding the calculation of p(x) for the 2-D space, we 
start with the 2-D free-space Green's function {k^r) 
|20j . For the TM polarization, the far-field electric field 
is given by 

E^(r)<x J As(r')E m (r')H^(k \r-r f \)d 2 r'. (A13) 

Using the asymptotic form of the Hankel function given 
by Eq. (A4), we have 

E^(r) oc -L f Ae(r')E m (r')e^ kof - r ' d 2 r'. (A14) 
V r J 

Comparing Eq. (A14) with Eq. (8), one can expect that 
for the 2-D case, the geometric integral in Eq. (33) would 
be (for the TM polarization) 



(A15) 



which is just Jo(fco|a:n — ^n'|)- 

A final note on the transition from Eq. (19) to Eq. (21). 
The simplification takes advantage of the following prop- 
erty of the spatial part of the electric field: 



V x (V x E{r)) - u 2 c ne(r)E{r) = 0. 



(A16) 



Obviously, E m (r) satisfies the above equation, but 
[Ej(r)] do not. In this paper, for [Ej(r)], we have ap- 
proximated e(r) as eq, i.e., the effect of the microres- 
onator to the Green's functions has been neglected, which 
is the essence of the volume current method that results 
in much simplified solutions with acceptable accuracies 
(a more physical explanation is that a proper choice of 
radiation modes should ensure that they are orthogo- 
nal to the WGM modes so they do not couple to the 
WGM modes without scatterers. However, [Ej(r)] given 
by Eq. (18) do not satisfy this property). Therefore, a 
more rigorous calculation should take the effect of the 
microresonator into account, and in Eq. (33) the exact 
Green's functions have to be used. For 2-D cases, such 
a task is relatively easy, and one illustrative example is 
provided in Appendix B. However, for 3-D cases, the ac- 
curate calculation of the Green's functions is usually dif- 
ficult 1201 ■ 
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Appendix B: 2-D FEM Simulation and Discussions 

In this part, we will describe the simulation method to 
obtain the complex eigenfrequencies of a 2-D microdisk 
resonator with scatterers attached on its surface. We 
use the finite element method (FEM), which is gener- 
ally much faster compared with the finite-difference time- 
domain (FDTD) method. To obtain the scattering loss, 
perfectly matched layers (PMLs) are implemented based 
on the stretched coordinate method [29 . For example, 
in the cylindrical coordinate system, for the TM polar- 
ization, we have [30 



The memory requirement is not that demanding (desk- 
top computers with 8 Gb memory run the simulation 
smoothly) , and each simulation typically takes a few min- 
utes. Moreover, the scan of the position of the scatterer 



d d_ 1 d 2 



E z 



-k 2 n 2 s 2 E z 



(Bl) 



where n is the refractive index at each region; s p is the 
complex coordinate stretching factor for the PML (light 
is only attenuated in the increasing p direction); and k 
is the eigenvalue we try to obtain, which is related to the 
complex eigenfrequency u as U = fcoc. 

Equation (Bl) is implemented and solved as a partial 
differential equation (PDE) in the commercial software 
COMSOL 26] . Because COMSOL does not provide the 
cylindrical coordinate system for structures without axial 
symmetry, Eq. (Bl) is converted back to the Cartesian 
coordinate system as 




dx 2 



dy 2 



E, 



= -kln 2 s 2 E 7 



(B2) 



FIG. 9. Simulated structure in COMSOL. The two small 
scatterers have been exaggerated in size for the illustration 
purpose. 



can be facilitated with the use of the COMSOL Script 
(or COMSOL with MATLAB) [2B] . 

As a supplementary discussion, here we provide a brief 
derivation for computing the accurate p{4>) in the 2-D 
space. Instead of using the approximate Green's function 



Figure 9 shows the structure we simulate, where the mi- #o ( k o r ) whi ch neglects the effect of the microresonator 



crodisk is centered at the origin. s p is chosen to be the 
following form for the PML region: 



(see discussions at the end of Appendix A), we seek the 
exact Green's function by solving the following equation 
(for TM-polarization)[2"0]: 



s p = 1 



d 2 



(B3) 



pd P yP dp> 



1 d 2 

p 2 d(f> 2 



+ kin 2 



where po and d are the starting radius and the thickness 
of the PML region, respectively, and a is a parameter 
that can be adjusted for the optimum performance of 
PML (we take a = 3 in our simulation). s p is 1 for other 
regions. 

Similarly, for the TE polarization, the equation can be 
implemented as 



g(r,<fi) = 6(r-R)5(<p). 

(B5) 



d_ 



d 

n 2 dx 



d 
dy 



d 
n 2 dy 



— k s p H z 



(B4) 



The placement of n 2 inside the first-order derivative is 
to ensure correct boundary conditions when the PDE is 
solved (i.e., E$ to be continuous) [3T] . 

For FEM simulations, mesh size has to be small enough 
to avoid artificial effects. In our case, cubic meshes with 
grid size less than 50 nm are employed (of course, for ar- 
eas surrounding the scatterers the mesh has to be finer). 



The solution has the following form [TT] : 

g{r, 0) = J2 a n Jn(k n d r)e m ^ for r < R, (B6) 

n 

g(r,<l>)=Y / bnHW{k n r)e in t for r > R, (B7) 

n 

where n d and n Q are the refractive indices of the mi- 
crodisk and the surrounding medium (air here), respec- 
tively, and a n and b n are the corresponding expansion 
coefficients. Applying the boundary conditions of g(r, 4>) 
at r — R gives 



a n J n (k n d R) = b n H\ >(kon R), 
k n d a n J' n (k n d R) - kon Q b n H'^ l \k Q n R) 



which yield 



2tt' 



(B8) 
(B9) 
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J n (kon d R) 



n d Hn \k n R)Jn{kondR) - n Hn 1 \k n R)J n (kondR) 



27rfcn 



(BIO) 



Inserting the obtained Green's function into Eq. (33), one 
can easily derive p(4>) for the TM-polarized light as 



(a) 



p{4) 



• cos no 



(Bll) 



where <j> is the angular separation between the two scat- 
terers as defined in the inset of Fig. 4. The numerical 
result of Eq. (Bll) for the example studied in Fig. 4 is 
shown by the blue dashed line in Fig. 6. 



Appendix C: Scatterer Distribution 

In this part, we will construct a distribution of small 
scatterers that satisfies the statistical rule given by 
Eq. (96). We consider identical scatterers, with the shape 
shown in Fig. 10(a). The height a and the width W 
of each scatterer are assumed to be much smaller than 
the wavelength. We divide the perimeter of the mi- 
crodisk resonator by N divisions to allow for N scatterers 
(W = 2irR/N), and each scatterer can be either point- 
ing outward or inward, with a parameter x k being +1 for 
the former and —1 for the latter for the fcth scatterer. A 
set of [x k ] (k = 1, 2, • • • , N) then describes the scatterer 
distribution on the perimeter of the microdisk resonator 
, and 

Ar((j>) — ax k , with k = Round (-=—), (CI) 

where RoundQ denotes the nearest integer function |25j . 
We generate [x^\ using the following statistical rule: 

x k e {-1, 1}, P(x k+1 = x k ) = \{l + X), (C2) 

where x f &us m the range between (0, 1). It is easy to 
verify 



E(x k x k+n ) = x lnl , Vfc, 



(C3) 



which states that the correlation between x k and x k+n 
only depends on their index difference. From Eq. (CI), 



< Ar((f>) Ar(0 + 0') >= a 2 E(x k x k+n ) sa a 2 X ^ , 

(C4) 

where we have approximated the index difference n cor- 
responding to qb' phase shift as n s» (f>'N/2ir, which is 
valid when |n| is much larger than 1. As a result, 



< Ar{4>)Ar((f> + <p') >w a 1 exp - In (x _1 ) 



2tt 



( b ) [\F(k n )\] ( c ) [Angle (F(k n W* 




0.2 0.4 0.6 
1^)1 



I 1 

Angle (F(*jo))/ir 



FIG. 10. (a) The scatterers considered in this model are all 
identical, and can point either outward (+1) or inward ( — 1) 
on the surface of the microresonator.(b)-(c) Amplitude and 
phase values of one simulation result of [F(k n )] (R— 10 /im, 
N = 5000, x = 0.926 so L c — 160 nm. For simplicity, we have 
assumed Sn 2 a = 1 in Eq. (C7)). Strong variations of [F(k„)] 
with the index n can be observed, (d)-(f) The value for a 
specific F(kn) (n — 90) is recorded for multiple generated 
[xk] (20,000 runs). The histograms show that the amplitude 
of [F(k n )] follows a Rayleigh distribution (Fig. 10(d)) and 
the phase of [f (fc n )] follows a uniform distribution in (—it, tt) 
(Fig. 10(f)). Figure 10(e) plots of the cumulative distribution 
function of the obtained |F(fc n )| (blue solid line) versus a fit 
for the Rayleigh distribution (black dashed line), where the y 
axis is shown using the logarithmic scale. 



Comparing Eq. (C5) to Eq. (96), we find they have sim- 
ilar expressions (remember x = R<j>), and the correlation 
length L c can be identified as 



2ttR 



W 



iVln^- 1 ) hi Or 1 ) 
[F(k n )] can then be calculated from Eq. (59) as 



(C6) 



F(k n ) 



2it5n 2 a 
N 



E 



x k e 



—il-rrkn/N 



R 



-X n , (C7) 



(C5) 



where [X n ] are the DFT (discrete Fourier transform) of 
[x k ]- 

Numerical experiments are performed in MATLAB by 
generating [x k ] based on Eq. (C2) and computing [F(k n )] 
from Eq. (C7). Figures 10(b)- (c) show one example of 
[F(k n )], which confirm that there are strong amplitude 
(Fig. 10(b)) and phase (Fig. 10(c)) variations with the 
index n. In Figs. 10(d)-(f), we monitor the value of one 
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F(k n ) (n — 90) for each generated [xk], and the his- is a Rayleigh distribution while the phase distribution of 
tograms imply that the amplitude distribution of F(k n ) F(k n ) is a uniform distribution in (— n, ir). 
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